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Abstract 


The one-dimensional Burgers model of turbulence is investigated by com- 
puting the functional integral expression for the correlation function* based on the 
Hopf theory of statistical hydromechanics, with the aid of a high-speed computer. 
The initial probability distribution of the velocity is assumed to be normal with 
zero mean and with a Gaussian covariance function. The manner in which the 
energy decay curve changes under variation of the Reynolds number R implies 
the existence of a certain asymptotic curve for R -* 00 . The values obtained for 
the correlation function at some instants indicate that the inverse-square law 
for the energy spectrum holds in some wave-number range for high values of R. 
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1. Introduction 


Our last report 1 sketched briefly the method for computing the correlation 
function of the turbulent velocity in the Burgers fluid model, 2 using a functional 
integral expression. Here a full treatment of this method is presented as well as 
further results. 

It is well known that the statistically considered Burgers equation offers 
insight into some of the properties of real turbulence. Our approach is to 
investigate the Burgers model on the basis of the Hopf theory of statistical 
hydromechanics, 3 with the aid of a high-speed computer where necessary. 
According to Hopf, the probability distribution of the velocity field u(x) at time t 
is given by the following characteristics functional: 

t) = exp 

where Uq and y are the scalar fields in x {- m < x < + a>); T t u 0 (x) is the field 
which at time t has developed from the field u Q (x) at time t = 0 according to the 
dynamical rule given by the Burgers equation; P 0 (u 0 ) denotes the probability 
distribution of u 0 at t = 0; and Jq 3 1^ denotes a Lebesgue integral with the 
measure 3^ over the whole set ft composed of all realizable functions u Q (x) . 
Formula (1.1) is meaningful as a solution to the turbulence problem, once P Q is 
specified at the initial time and provided that the right-hand side is calculable, 



y (x) T fc u Q (x) dx> SP Q (u Q ) , 


( 1 . 1 ) 
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since (1.1) satisfies the Hopf equation for characteristic functionals and the 
Hopf equation is completely equivalent to the well-known hierarchy of coupled 
equations for the various-order simultaneous correlation functions of the turbu- 
lent velocity field. It is easy to obtain the formula for a correlation function 
from (1.1). For example, 


S 2 

i 2 Sy^) Sy(x 2 ) 


y = 0 


= Tt u 0 (*i) TX (x 2 ) SP 0 (li ( 


= <u (Xj) u (X 2 )> t , 


( 1 . 2 ) 


which is the two-point correlation function at time t. Here 6 <b/8 y(x) denotes 
the functional derivative 3 - 4 of 0 at x. 

The key to making the right-hand side of (1.1) or (1.2) calculable seems to 
be to use a properly defined functional integration 5 to replace the abstract 
Lebesgue integral over a function space. Let us formally set 


sP o( u o) - P ( U 0> °) Su o- U- 3 ) 

Here p can be understood as the probability density for the function u 0 (x), when 
u Q (x) is expanded in terms of the real-valued orthonormal function set { s k (x) } 
according to 

2N 

u f O) = a k s k< x )’ (!- 4 ) 

k a 1 

and 8 u 0 is read as 


4 


» 


(1.5) 



2N 


k - 1 


da k 


provided, of course, that 



( 1 . 6 ) 


The larger the value of 2N, the more variety our function space can include. 

The most ideal situation would thus be reached in the limit as 2N-* This limit- 
ing operation cannot be carried out in (1.5), because the infinite-dimensional 
volume element is meaningless, 6 but should be applied either to (1.3) or to the 
integral in (1.1), As long as this is kept in mind, we may safely use the symbolic 
notations which are obtained by taking the superscript 2N away from both sides 
of (1.4) and (1.5). Practically, however, we attempt to approximate (1.1) or (1.2) 
with 2N < , as described in the next section. The accuracy of approximation 

therefore depends on the rate of convergence to the limit, which may in turn 
depend on what orthonormal function set is used. 

Thus, if T^ (x) is known explicitly in terms of u 0 and if p(u Q , t 0 ) is properly 
assumed, we find no essential difficulty in proceeding with the calculation in- 
volved in (1.1) or (1.2) except that a high-speed computer is generally necessary 
to estimate the functional integrals needed. Indeed, for the Burgers model, Hopf 7 
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and Cole’s 3 general solution can be used as T t u 0 (x). In the next section, the 
method of calculation is developed in detail on the premise that the initial probabil- 
ity distribution of the velocity is normal. This premise leads to the following 
form of initial characteristic functional: 


0(y„» °) = ex P y- f u ( x )y 0 ( x ) dx - \ f f Q( x > x Oy 0 ( x )yo( x< ) dxdx/ ) 


(1.7) 

where it is easy to know that the functions TJ and Q are the average and covariance, 
respectively, of the stochastic velocity field Uq (Cf. (1.2)). We assume further 

U (x) = 0 and Q(x,x') = e -(x - x '> . (1.8) 


i.e. the average velocity vanishes and the initial correlation function is spatially 
homogeneous and Gaussian. In order for the condition (1.7) to be consistent with 
(1.1) at t = 0, the function p(u 0 , 0) defined in (1.3) must be 


P (u 0 , 0) = 



u 0 ( x )y 0 00 dx 


0(y o * °) § y 0 - 


(1.9) 


(y n and S y are understood in the same sen$a as described for u and SuJ, 

'* 7 0 J o o o 7 

The assumption of normality of the initial probability distribution is made because 
it can be proved from information theory 9 that this distribution is the most unbiased 
one, i.e. it makes the information entropy a maximum, if we have no certain 
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information other than J and Q. (See Appendix.) Some departure from normality 
may conveniently be considered by introducing a functional Gram-Charlier series, 
and more general assumptions for U and Q may also be used if necessary. How- 
ever, the present research is confined within the limits set by the above premise, 
since our main purpose is to verify the practicability of our fundamental approach. 


2. Method of Calculation 


i) Diagonalization of the normal probability density 

.... . ■■■ — i —• — ’ — I..' i i ■ — ■ ■ ■ ■■W— ■■ ■ « '■ 

We first calculate (1.9) in order to obtain a simple explicit form of (1.3). 
This can be achieved by a transformation of the coordinate frame of the function 
space. To correspond with (1.4), let us start from the expression 


yf (x) = 


L 

i = 1 


a. s. (x) 


( 2 . 1 ) 


At x = x (m = 1, . . . , 2N), we set 



N 


Y~' i ^( k j ) c °sk.x m + i (k.) sin kjX m > 
j = l 


Ak. 

j 


( 2 . 2 ) 


where A k. is the length of the j th interval when the number range (0, <x>) is 
divided to N intervals , and k. is the representative value of k in thej th interval. 
The 2N simultaneous equations obtained by equating the right- hand-sides of (2.1) 
and (2.2) for each value of m give a one-to-one correspondence between (a t , . . . , 

a 2N ) and (T? < k i ) ( Ak i) 1/2 » •'. • r^(k N ) k N ) 1/2 » £ (kj) (Akj) 172 , . . . , i (1%) 

(Ak^) 1/2 ). If we multiply these 2N equations on both sides by s n (x m )A x rii , where 
A x m is an interval about the representative point x m of such length that the sum 
of all disjoint A x m covers the whole x space, and if we then sum over m, we 


obtain 


cos k. x dx 
) 


:.) (* s 

j J | n 

J- CO 


+ £ (k.) s (x)sinl xdx> + e 


(2.3) 


by the properties of orthonormal functions. Here e is a quantity which tends to 
zero as all Ax - 0 (i.e. as 2N~ ). The Jacobian of the transformation is 

calculated as follows: 


J^ck.KAk.)" 2 } 


Ak.\ 1/2 r 00 

— -j I s.(x)cos k.xdx+0(d) = t ijl (2.4) 


-dU(kj)(Ak.) 1/2 } 


MkA 1/2 r 

= i s. (x) sin k. x dx + 0 (e) = t. . . 

W / L J 1,N + J (2.5) 


Hence we have 


Y\ (‘ij ‘-ti tti.N.j U,N + i ) =7 dxdx' S, (x) S{(x') 

j _ j J— CD J—OD 

N 

• ^ (cos k^x cos k^.x' + sin kjX sin kjX') Akj + C (e) 

j = l 


= ^ + 0(e) + e' 


( 2 . 6 ) 
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where is the Kronecker symbol and the quantity e' tends to zero as all 
A k. -* 0 (i.e. as N °°). This indicates that the tr ansformation of the represen- 
tation of the function space from (a t , . . . , a 2N ) to ( 77 ^ ) (Akj ) l/2 f . . ., 7 ? (1^) 
(Ai^)^ 2 , £ (k^ ) (A k l ) i/2 , . . . , £ (Ak^) 1 /2 ) is orthogonal in the limit N -> “ . 

Accordingly, 8 y 0 may be expressed in our symbolic sense as 

5 y 0 = ~jj~ d7 7( k i> 


’’dkA l ' 2 /dk. \ 1/2 

.17) ‘Wfe) 


(2.7) 


At the same time, (2.2) reduces for almost all x to 


y 0 ( x ) 



cos kx + £ (k) sin kx} dk . 


( 2 . 8 ) 


If we introduce the relations 


77 (k) 


_ z (k) + z (-k) 
2 1/2 


i (k) = 1 < z ( k ) - 2 (~ k )> 
2 i/2 


(2.9) 


(2.8) may be rewritten as 


-to 

y 0 (*> - z (k) e ikx dk , 

(277) 1/2 I* 


( 2 . 10 ) 


where z*(k) = z (— k) , the asterisk denoting the complex conjugate. With the aid of 
these relations the quadratic form in (1.7) can be rewritten in a diagonal form 
as 
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1 

2 


cr (k) z (k) z (-k) dk 



Q(x -x') y Q (x) y 0 (x') dx dx' 




1 

2 



Re cr (k) {t] 2 (k) + £ 2 (k)} dk 


(2. XI) 


where Re indicates the real part and cr (k ) is the Fourier component of Q(x), 
which when calculated from (1.8) is found to be 


cr(k) 


£ 


QOO 


: ikx dx = 7 r l/2 e' k2/4 


( 2 . 12 ) 


To make possible the calculation of the functional integral in (1.9), we introduce 
a modified Fourier expansion of u 0 (x) as follows: 


Uo(x)= ^I 


[Re cr (k)] 1//2 {A (k) cos kx + M (k) sin kx } dk . . (2.13) 


Hence, we can compute 


Vo ( x ) u o ( x > dx = [Re cr (k) ] 1/2 .{A (k) 17 (k) + M (k) £ (k) } dk , (2.14) 

J— CO Jo 


and finally can arrive at the formula 


p (u Q , 0 ) = exp 


-7 21 {A2(k i )+M2(k i )>dk i Tf (r^oT)} ■ < 2 - 15 > 
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Now 5 u 0 may be expressed in the same way that S y 0 is expressed in (2.7). If (2.13) 
is compared with (2.8), it is seen that the proper expression is 





dA(k j ) 


dk. \ 1/2 

~ ) [Recr(k.)l 1/2 dM (k.) 



1/2 


(2.16) 

These last two equations now provide a simple, explicit, direct-product form for 
(1.3): 



(2.17) 


which demonstrates that the stochastic variables { A(k. ) dkj /2 , M(k .) dkj /2 } 
distribute standard-normally, independently of each other. It is worth noting that 

(2.17) is just the Gaussian measure on the function set of A(k) cross M(k) 
established by Friedrichs, 10 so that (1.1) or (1.2) is just a functional integral 
with the Gaussian measure, the integrand of which is a functional of A (k) and 
M(k) through the relation (2.13). 
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ii) Monte Carlo quadrature 


For the Burgers fluid model, T* u Q in (1.1) or (1.2) can be written as 



(2.18) 


(the, Hopf- Cole solution 7, 8 ); this satisfies the dimensionless Burgers equation, 


3u _ B u 1 B 2 u 
Bt U Bx + R 3x 2 


(2.19) 


where R is the Reynolds number referred to the initial characteristic correlation 
length and the root-mean-square value of the initial velocity (Cf. (1.8)). 


In order to estimate the functional integral (1.2) for any values of x t , x 2 , and 
t, it is most convenient to use the Monte Carlo quadrature for a multiple integral 11 
on an approximate basis such that A k. and N are kept finite in a proper way. 

That is, we will use the so-called cylinder fu ctional approach. Estimation of 
(1.2) then reduces to the averaging of many sample values of T t u 0 (x x ) T* u Q (x 2 ) 
(called the estimator 11 ) which corresponds to taking many sets of 2N independent 
standard-normal random numbers as values for the variables {A (k.)dk 1/2 , M(k. ) 
dkj /2 > . 
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Thus, (2.13) may be replaced by 


u 0 (x) = 


77 


1 

1/4 


IN 

L 

j=i 


n 

-k: / 8 


{A(k.) cos x 


+ M(k.) sin k.x} Ak; (2.20) 

here (2.12) was substituted for & (k) and all A k. were put equal to A k so that k . = 

(j - 1/2) Ak . The maximum wave-number k c to be considered for u 0 (x) is 
given by k c = NAk. The value of k c , which is chosen so as to make the error of 
(2.20) due to the neglect of higher wave components less than 1%, may be determined 
from the equation 


{-Hi 


e" k /8 dk 


= 1 / 100 , 


( 2 . 21 ) 


which yields k c = 5.17. Once k c is fixed, Ak depends on N. The question of 
how large N should be was discussed in our last report; 1 results presented herein 
are given for values of N as large as 30 and even 45 to insure convergence. The 
larger N = k c /Ak is, of course, the smaller is the error due to using the rec- 
tangle rule as an approximation to the integral expression (2.13). Using (2.20), 
the integral over x" in (2.18) may be expressed as 


r 

J 0 


Uo( x") dx" =_L 


N - k. / 8 


L 


k. 


(A(k. ) s in kj x' 


- M(kp cos kj x' +M(kp} Ak 
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% 

But, in the case where Ak is not necessarily small, the alternative expression 



u 0 (x") dx" 




(kj ) s in kj x' 


-M(k.)cosk. x'] sin Ak x'/2 + Mrk ^ Ak 
J 3 J x'/2 3 


( 2 . 22 ) 


is much more efficient. Indeed, the factor sinA k (x f /2)/(x*/2) prevents the 
deterioration of the Monte Carlo quadrature caused by a large vibration of sample 
values for a large value of x\ The integrals of both the numerator and the de- 
nominator of (2.18) were calculated by the Simpson rule, using (2.22). In this 
calculation there are two points to be noted. First, the exponential function in 
the integrands should be normalized to a value as close to its maximum value 
as is possible, because at some times it becomes enormously large and at other 

4 

times infinitesimally small. Second, the interval of integration cannot be infinite 
in practice, so that integration is stopped at the points where the integrands 
decrease to less than 10' 7 times their maximum value, taking account of the 
rather rapid, monotonic decrease of the factor exp { -R(x - x’) 2 /4t } for [ x r | 
large. For convenience, a new variable, s = (R/4t) 1/2 (x - x’), was introduced 
in place of x*. The accuracy of our Simpson numerical integration may be 
judged by the fact that tb~ value of the integral was unchanged except for the 
fourth digit when the integration step length was decreased by a factor of 2. 
Examples of the calculation of (2.18) are given in Figs. 8 and 9. 


15 



Our standard-normal random numbers were generated by the simple but 
highly reliable method of Kronmal. 12 The number of sample values for the 
estimator in our Monte Carlo quadrature was between 300 and 500. The 
reliability of this quadrature may be judged from the examples shown in Figs. 
2 and 4. 
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3. Results and Discussion 


i) Features of the energy decay 

In Fig. 1, twice the calculated dimensionless energy of turbulence, <^[u(0)] 2 ^ , 
is plotted against the dimensionless time (normalized with the initial characteris- 
tic correlation length divided by the root-mean-square of the initial velocity field) 
for various Reynolds numbers. The Monte Carlo quadrature was done with N = 

30, and the number of samples is 300. The solid lines show the corresponding 
results from the linear theory, <([u(0)] 2 )> t = 1/(1 + 8t/R) l/2 , where the non-linear 
term in the Burgers equation is neglected. The dotted line gives the correspond- 
ing results from the approximate theory of Meecham and Siegal 13 for R = 100. 

The reliability (or stability) of our Monte Carlo quadrature is indicated by Fig. 2, 
which shows how the value of the integral depends on the number of samples taken 
in the case R = 100. 

Considering the approximate nature of Monte Carlo quadrature, the agreement 
between our results and the linear theory for R = 0.1 and 1 is excellent at the 
initial time; but after the energy has decayed to one-tenth of its initial value, our 
results begin to deviate from those of the linear theory, giving rise to a little 
more rapid decay. This seems to be due to a nonlinear effect which transfers 
energy from the small wave-number components to larger wave-number com- 
ponents, where viscous dissipation is stronger. 

For R greater than 10, the character of the energy decay is quite different 
from that expected on the basis of the linear theory. Even more interesting. 
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the plotted points suggest the existence of an asymptotic curve for R - a> . Indeed, 
the points plotted for R = 100 seem already to represent the features of the 
asymptotic curve to good accuracy. (Cf. those for R = 1000.) As can be seen, 
the theoretical curve of Meecham and Siegal deviates considerably from this 
expected asymptotic curve. Their approach is based on the Wiener-Hermite 
expansion of the random field u(x) , but with some drastic assumptions in order 
to avoid an open hierarchy of unknown functions and other complexities involved 
with the calculation, whose validity was questioned by Orszag and Bissonnette. 14 

In Fig. 3, the manner in which the result depends on N is shown for the case 
of R = 100. The data for N = 15 were also presented in our last report; there, 
on the basis of a comparison with the data for N = 10, it was argued that the N = 

15 data accurately represent the solution in terms of the functional integral (1.2). 
This comparison, however, was not satisfactory because the values of k c for the 
two cases were considerably different. As can be seen in Fig. 3, the data for 
N = 15 are still rather far from the expected limiting values for N - °o with the fixed value 
of k c given by (2.21). It appears that the data for N = 30 are much nearer the 
limiting values. This indicates that the accuracy of the present calculation is 
higher than before. 

There are some arguments to the effect that the energy of the (spatially 
periodic) Burgers’ model turbulence should decay like t -2 in the limit R -» , 
if t is large enough. 15 Our results for large N include no symptom of decay like 
t“ 2 , at least up to t = 100. However, we do see symptoms of such a rapid decay 
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for small values of N. (See Fig. 3.) The smallness of N means that rhe initial 
velocity fields are restricted to a considerably small class of the function space 
Q, as is known from (2.20). For the ensemble of such a relatively simple class 
of initial velocity fields, the saw-tooth wave argument 15 appears to be easily 
realizable, because (2.20) ';hen has the nature of a velocity field which is strongly 
spatially periodic (relative to the initial characteristic length of correlation). 

But the larger N is, the larger the class of initial velocity fields involved is, and 
the more difficult it is to realize the so-called saw-tooth wave argument, in a 
short time, at least. Here it is worth recalling that there exists a similarity 
solution of the Burgers equation which has the shape of a solitary wedge with the 
hinge point fixed on the x axis, and that the energy of this velocity field decays 
not like t“ 2 but rather like t~ 1/2 . 2 If tbi 6 type of velocity field comes into the 
ensemble, obviously the rapid decay based on the saw-tooth wave argument should 
be reconsidered. There might be various other solutions with different decay 
patterns. In the limit as N - ro the ensemble should include all these solutions. 
This fact seems to give theoretical support to the present result of a rather slow 
decay of the energy. It is most interesting to note that in the later period of 
decay (t > 3) our result is close to Burgers’ prediction, 2 based on a statistical 
assumption, of a t _2/3 dependence. Burgers assumes that during the intermediate 
period of the decay, the velocity field over the infinite domain may be approxi- 
mated for large B by a series of discontinuous straight-line segments of positive 
slope which decreases with time but of random length and random magnitude 
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of drop from one segment to another. He also takes into account frequent 
coalescence of discontinuous points, i.e. shocks. 

The final result obtained of course would depend on what ensemble of 
velocity fields is taken initially, i.e. what is prescribed as <£(y 0 , 0). Our choice 
is indicated by (1.7)— (1.8) ; and, as we have already said on an information- 
theoretical basis, this ensemble is the one most often realizable, unless we have 
other peculiar information. If N is too small, it is obvious tU at the calculation 
cannot give an exact result for this initial ensemble, since in that case the 
functional integration has been carried out imperfectly. 

ii) Correlation function and energy spectrum 

The correlation function <(u(0)u(x)^ was calculated at t - 1 for R = 1000, at 

t = 1 and 3 for R = 100 and at t = 1 for R = 1, as is shown in Fig. 4. The solid 

2 

line is the curve at t = 0, i.e. Q(x) = e~ x . The dotted lines are the theoretical 
results of Meecham and Siegel 13 at t - 1 and 3 for R = 100. The reliability of 
our Monte Carlo quadrature is observed in Fig. 5, where the number of samples 
was taken up to 500. 

• Fig. 4 shows the general trend of the correlation function to become flatter 
and flatter as time goes on. As expected from the observation in the preceding 
subsection, there is so little difference between the correlation functions at the 
same instant t = 1 for R = 100 and R = 1000 that they are considered to represent 
the asymptotic curve of the correlation function for R -* co . However, the 
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correlation function at the same t for R - 1 is very different from that for a 
higher Reynolds number; it becomes flat faster than the latter, because of the 
stronger viscous dissipation. It may be seen from the figure that the correlation 
function at t = 1 for R = 1 is quite similar to that at t = 3 for R = 100, which con- 
tains a comparable amount of turbulence energy. But the fact that these two 
curves cross each oth</r somewhere between x = 1.0 and x ~ 2.0 with different 
curvatures seems to suggest some delicate structural difference between them. 
This may be clarified by comparing their energy spectrums. 

The energy spectrum corresponding to each correlation function can be 
calculated from the equation 


F,(k , t) = 



cos kr *Q(r, 


t) dr, 


(3.1) 


where Q(r, t) = <(u(0)u(r))> t . The curve of Q(r, t) was made by connecting the 
plotted points in as smooth a way as possible, but the result of numerical calcu- 
lation for E(k, t) (using Filon's method 16 ) predicted not a smooth curve but a 
somewhat rugged curve, as shown in Fig. 6. This seems to be an unavoidable 
effect caused by the error inherent in Monte Carlo quadrature. Nevertheless, 
it is interesting to note that there is an appreciable region of the energy spectrum 
in which the slope is very nearly k~ 2 for R = 100 at both t = 1 and t = 3. It is 
evident from the figure that this well-discussed feature of the energy spectrum 
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4k 

of the Burgers' model turbulence 2 ’ 13 ’ 1S ’ 17 appears in a very short time and is 
kept pretty long for R = 100, while the same feature is barely observable 
at t = 1 for R = 1. This contrast between the energy spectrum at t = 3 for R = 100 
and that at t = 1 for R = 1 explains to some degree the above-described structural 
difference between the correlation functions for these two cases. We may under- 
stand this fac L from the nonlinear effect of the Burgers equation; indeed, without 

the strong nonlinear interaction the energy in high wave-number regions would 

2 

decay almost according to e _2k t/R , while if the energy transfer from lower to 
higher wave-number regions of the spectrum occurs through the nonlinear inter- 
action, the efiect of the relatively strong energy decay of the high wave-number 
components will be considerably mitigated so that the k~ 2 spectrum law may be 
recovered in some range of k, which may be called the inertial rmge. 

As shown in Fig. 4, Meecham and Siegel's theoretical results are different 
from our results at t = 1 and 3 for R = 100; but from the viewpoint of the energy 
spectrum, shown in Fig. 6, both have a somewhat similar slope, except that 
their theory predicts a smaller energy of turbulence at t = 1. The same tendency 
is also seen in the numerical experiment performed by Jeng et al, 18 although their 
result is not plotted in the figure in order to avoid complication. It may be noted, 
however, that their experiment is very similar to the present calculation, in 
that both deal with the operator T* in (1.2) exactly and both try to compute 
an average property for an ensemble of initial velocity fields by sampling. But 
Jeng et al. limited the ensemble to one special type , made up from n independent 
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uniform random numbers (n seems to be 50, judging from an illustration in their 
paper), so that ef> (y 0 , 0) in their case has a form quite different from that of (1.7). 
Furthermore, in order to calculate the statistical properties of the velocity field, 
they appealed to the operation of space averaging (assuming statistical homogenity 
and the so-called ergodic hypothesis) rather than averaging over a large number 
of members of the ensemble. These facts give their experiment a flavor quite 
different from our approach. If they had simply proceeded with the latter 
averaging operation, they would not have needed any hypothesis beforehand 
but, after making sure of the reliability of their samplings, s hoi, Id have been 
able to prove it. In spite of these rather great differences between their method 
and ours, it is notable that the k~ 2 spectrum law is qualitatively recovered by 
both. It would be even more interesting if it were known what result they had 
for the feature of energy decay, because a comparison with Fig. 1 would then 
show how the result is affected by the initial condition on 0 . 

It is worth noting that all energy spectrums in Fig. 6 tend to almost the 
same value as k 0, regardless of t and R. This value may be estimated from 
the theoretical curve at t = 0, giving 


E(0, 0) =i 

77 - 



dr = 77 1/2 /2 . 


(3.2) 


This fact is consistent with the existence of the so-called Loitsianski constant 
for Burgers’ model of turbulence. 2, 13 We have verified that at k - 0 our results 
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for different t and R recover the value (3.2) with a relative error of less than 
5%. This may be considered as another measure of the order of accuracy of our 
Monte Carlo quadrature. 

iii) Deviation from normality 

It is hard to make a decisive quantitative statement on the extent to which 
our probability distribution functional of the velocity field deviates from normality, 
because we must always take into account the error intrinsic to Monte Carlo 
quadrature. To determine this deviation, we must deal with a correlation whose 
order is higher than two, and the higher the order of correlation dealt with, the 
larger would the value of N be needed in order to keep enough accuracy. Fig. 7 
shows the result of calculating, at N = 30, the curtosis for many values of R and t. 
In our case, the curtosis at t = 0 should be exactly 3, since we have started from 
a normal distribution. The only plausible conclusion, therefore, is that there is a 
tendency for the curtosis to decrease very slightly with increasing time, i.e. for 
the distribution to become slightly flatter, at least for R greater than or equal to 
1. Calculated values of the skewness factor are distributed almost randomly about 
zero; hence it may be concluded that there is no significant deviation from the 
zero value. 

iv) Time-development of an individual velocity field 

Finally it is worthwhile to look into the behavior of a typical individual 
velocity field. Fig. 8 shows the time development, for R = 100 and N = 30, of the 
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particular velocity field which at the initial time was made up from a certain set 
of standard normal random numbers, plotted over the interval [-2.5, 2.5] . 
(Actually, the random numbers in this case are taken between the 61st and the 
120th of those generated by Kronmal's method. 12 The general feature of 
an individual velocity field was unchanged even if another set of random numbers 
was taken.) At t = 1 a shock-like sharp discontinuity can already be observed; it 
moves along according to. its own dynamics, 2 but diminishes gradually because 
of the relatively predominating (with decay) viscous effect. What Burgers 2 calls 
the hinge point around which the curve of u(x) rotates in a clockwise fashion with 
time, can be seen at x = 2.3. Another hinge point at x = -3.6 is predicted. 

The slope of u(x) appears to be almost (t + t Q ) _1 near the hinge point and, for t > 1 , 
to be so everywhere except near a rounded shock-like region. (In this 
example t 0 = 1.) 

For high Reynolds number, these features identify the reality of Burgers' 
conception of the development of an individual velocity field pretty well. The 
period with t > 1 seems to be suitable for application of his simplified statistical 
assumption except for the rounded shock-like regions. For comparison, the 
calculated time development for R = 1 of the same initial velocity field is 
shown in Fig. 9. In this case the diffusive action due to viscosity predominates 
for all periods of time. The contrast with the former case is very remarkable. 
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4. Conclusion 


This work verified that the functional integral expression for the correlation 
function of the turbulent velocity field, derived from the Ilopf formalism for 
turbulence, 3 is directly calculable through the introduction of the functional 
integration technique given by (1.3)-(1.6) and the use of Monte Carlo quadrature. 

The present numerical result for the Burgers' model turbulence in the infinite 
domain of x is accurate enough to reveal a reasonable transition in the features 
of energy decay as R changes. For R less than 1 the turbulence energy decays 
just as the linear theory predicts over a considerable period of time, while 
for R greater than 1, the decay pattern changes rapidly with R so that even for 
R = 10 it is almost in accord with the predicted asymptotic (R - oo) pattern (at 
least until t = 100). Although the time-dependence of the asymptotic decay 
pattern is not in accord with the results of other approximate theories, except 
for Burgers', the energy spectrum at some instants in some wave-number 
regions has a tendency to behave at least qualitatively according to the k -2 
spectrum law, which is well predicted by all these theories including Burgers'. 

When the operator T t is not known explicitly, as is the case with the 
Navier-Stokes turbulence, we seem to face some difficulty in extending this 
functional approach. But it can still be used if T l u is developed as 

t) = ^ b k (t) s k (x) 
k 
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(Cf. (2.1)), where both x and u are vectors, provided that the time-development 
of { (t) } can be computed. This means that the T* is given numerically. 

Thus it may be expected that the present method of computation on the exact 
ensemble-mechanical basis will further contribute to the checking of the validity 
of approximate theories and will give typical data for engineering purposes. 



Appendix 


The average velocity field U(x) and the covariance function Q(x, x‘) are 
defined by the following equations: 

U(x) = I u(x)p(u) Su, (Al) 

Q(x, x') = I (U(x) - u(x)} {U(x') - u(x')> p(u) Su, (A2) 

h 

where p(u) is the probability density for u in the symbolic sense. We also have 

1 = I p(u) Su . (A3) 

J n 

Under the conditions (A1)-(A3), let us try to make the information entropy a 
maximum, where the information entropy is defined by the limit of 

S 2N = - p(u 2N ) log p(u 2N ) Su 2n / 2N as 2N - co , (A4) 

h 

according to Shannon. 19 The Lagrange multiplier method for the variational 
problem leads to the equation 
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0 


= - Su (log p(u) + 1} + I a(x) U(x) dx 
+ J| fa (x , x') (U(x) - u(x)} (U(x') - u(x')} dx dx' - 


(A 5) 


Here a(x), b(x, x’) and c are the undetermined multipliers. As a result, 


p(u) = e c 1 • exp 


- f«, V 


[-J 


a(x) U(x) dx 


) (U(x) - u(x)} (U(x') - u(x' » dx dx' 


(A6) 


where b(x, x f ) is assumed to be positive definite. The Fourier transformation 
of this is, in general, 


IJ 


-I 


<£(y) = c ex P I A(x) y(x) dx - II B(x, x') y(x) y(x') dx dx' 


(A7) 


where A(x), B(x, x’) and C are related to a(x), b(x, x f ) and c, and should be de- 
termined in such a way that and p are consistent with the conditions (A1)-(A3). 
Remembering the properties of the characteristic functional 4 > , we can obtain 


A(x) = i U(x) 


A 


B(x, x')=(l/2)Q(x, x') > 


(A 8) 


C = 1. 


J 
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Figure Captions 


Fig. 1. Features of the energy decay 

Fig. 2. Average values for <([u(0)] 2 ^ t versus the number of samples at various 
instants, with R = 100, N = 30. 

Fig, 3. Dependence of the energy decay on N, with R = 100. 

Fig. 4. Correlation function <\i(0) u(x)^> t for various values of R and t, with 
N = 30. 

Fig. 5. Average values for <^u(0) u(x)^> at t = 1 versus the number of samples 
at various distances, with R = 100, N = 30. 

Fig. 6. Energy spectrums E(k, t) for various values of R and t, with N = 30. 

Fig. 7. Curtosis as a function of time for various R. 

Fig. 8. Time development of b velocity field for R = 100. 

Fig. 9. Time deve ; : -ent %, partiecja?. velocity field for R = 1. 
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